home *** CD-ROM | disk | FTP | other *** search
/ Internet Info 1994 March / Internet Info CD-ROM (Walnut Creek) (March 1994).iso / networking / applic / ntp / dart.txt < prev    next >
Text File  |  1991-09-29  |  48KB  |  1,043 lines

  1.      On Synchronizing DARTnet Clocks to the Millisecond - a Report
  2.  
  3.                              David L. Mills
  4.                          University of Delaware
  5.                             3 February 1991
  6.  
  7. 1. Introduction
  8.  
  9. This is a status report on the hunt toward synchronizing DARTnet clocks
  10. reliably to the millisecond. Included in this report are the results of
  11. a preliminary investigation into the timekeeping quality expected for
  12. DARTnet experiments. The investigation suggests that DARTnet routers can
  13. attain and sustain accuracies to the order of a millisecond and
  14. stabilities to the order of a millisecond per day. It also points out
  15. certain routing problems that can seriously undermine the confidence
  16. that this level of performance can be guaranteed for all experiments.
  17. This should be a high priority problem to resolve. Also in this report
  18. is a description of the fuzzball clock statistics mechanism together
  19. with file formats and data interpretation which can be used to calibrate
  20. experimental results after the fact in order to realize the ultimate
  21. timekeeping accuracy.
  22.  
  23. All operating DARTnet routers are now running NTP Version 2 with each
  24. other and with two NTP stratum-1 (primary) time servers, wwvb.isi.edu at
  25. ISI and dcn5.udel.edu at UDel. The fuzzballs are running NTP Version 3
  26. as described in the latest revision of the protocol specification found
  27. in the PostScript file pub/ntp/ntpv3.ps.Z on louie.udel.edu. While
  28. Version 3 is backwards-compatible with Version 2, some minor differences
  29. in status reporting exist. The result of this is that the Unix program
  30. ntpq, ordinarily used to probe NTP time servers for status, operates
  31. slightly differently with the Unix xntpd and fuzzball servers. (The main
  32. difference is that the tattletale indicator reported just after the peer
  33. ID is not always meaningful for the fuzzball.)
  34.  
  35. The investigation also happened to reveal certain design problems with
  36. the fuzzball time servers used by DARTnet to synchronize local clocks to
  37. Coordinated Universal Time (UTC). These problems were due to the
  38. increasing demands placed upon the primary servers by a rapidly
  39. escalating cadre of users and involved relatively high timing jitter due
  40. to interrupt latencies, queueing delays and so forth. These problems
  41. have largely been alleviated by subtle and intricate jinks in the code;
  42. however, the problem remains that time service is suddenly getting very
  43. popular. To that end I have incorporated access controls in the UDel
  44. time servers designed to help hold down incidences of casual abuse, such
  45. as many users from the same network ganging up on the same server. I
  46. intend to update ISI and the remaining fuzzball time servers shortly. A
  47. message to that effect has been sent to the NTP interest list along with
  48. guidelines that should help reduce friendly fire casualties in future.
  49.  
  50. 2. Retrospective Timekeeping Data
  51.  
  52. I have configured a dedicated fuzzball dcn2.udel.edu (actually, one that
  53. is to be shipped to MIT when the NEARnet NSS is installed there). This
  54. machine is internally synchronized to the UDel master clocks and
  55. continuously watches every interface on every DARTnet router, as well as
  56. selected fuzzball time servers. The data are reported to a file which
  57. can be retrieved via FTP. For most purposes this statistical luxury can
  58. be ignored, since I intend to watch the clocks as part of my regular
  59. snoop rounds. However, for those experiments where timestamps reliable
  60. to the millisecond are required, experimenters may need to calibrate
  61. extracted timestamps using retrospective offset measurements available
  62. in this file.
  63.  
  64. The binary data file is called stats.dat and is accessible via login
  65. tcp-test, password cerf (in the finest historical tradition). It
  66. consists of records of eight 16-bit words in little-endian, binary
  67. representation.  The records are in the format:
  68.  
  69.      +------+------+------+------+------+------+------+------+
  70.      | Day  |     ToD     | Code |    Offset   |Delay |Disper|
  71.      +------+------+------+------+------+------+------+------+
  72.    word 0      1      2      3      4      5      6      7
  73.  
  74. Day       NTP day, where day 0 is 1 January 1900; to obtain Modified
  75.           Julian Day (MJD), add 41317 to this number
  76.  
  77. ToD       time of day relative to UTC midnight (milliseconds), high-
  78.           order 32 bits followed by low-order 32 bits
  79.  
  80. Code      two-octet code, coded as follows
  81.  
  82.           0-7  high-order octet of the peer status word described in
  83.                Appendix B of RFC-1119 and ntpv3.ps, coded as follows
  84.  
  85.                0    configured entry
  86.                1    peer authentication enabled
  87.                2    peer correctly authenticated (truechimer)
  88.                3    peer reachable
  89.                5-7  peer selection status, coded as follows (the
  90.                     character shown following the code is the tattletale
  91.                     indicator shown by the fuzzball monitoring program;
  92.                     however, the ntpq monitoring program may tattle
  93.                     different characters)
  94.  
  95.                     0    rejected
  96.                     1 x  passed sanity checks (tests 1 through 8 in
  97.                          Section 3.4.3)
  98.                     2 .  passed correctness checks (intersection
  99.                          algorithm in Section 4.2.1)
  100.                     3 -  passed candidate checks (only the first 10
  101.                          peers on the candidate list are considered for
  102.                          clock selection in the fuzzball implementation)
  103.                     4 +  passed outlyer checks (clustering algorithm in
  104.                          Section 4.2.2)
  105.                     5 #  current synchronization source; max distance
  106.                          exceeded (in the fuzzball implementation this
  107.                          means the synchronization distance exceeds 1000
  108.                          ms, normally considered real raunchy time)
  109.                     6 *  current synchronization source; max distance
  110.                          okay
  111.                     7    reserved
  112.  
  113.                     Note that these codes are cumulative; a status of n
  114.                     implies passing all tests at level n, together with
  115.                     all tests at levels less than n
  116.  
  117.           8-15 peer ID
  118.  
  119. Offset    signed offset of peer clock relative to local clock
  120.           (milliseconds), high-order 32 bits followed by low-order 32
  121.           bits
  122.  
  123. Delay     signed roundtrip delay via peer clock to the root of the
  124.           synchronization subnet (milliseconds); note: under some
  125.           conditions this number can be negative
  126.  
  127. Disper    unsigned dispersion  via peer clock to the root of the
  128.           synchronization subnet (milliseconds)
  129.  
  130. Note that the file is fixed size, normally 512K bytes, and is zero-
  131. padded at the end. Ordinarily, this provides for about a month of data
  132. recording, following which I rename the old file and create a new one.
  133. At present, there is no facility for automatically creating a new file
  134. when the old one fills up. Trust me.
  135.  
  136. 3. How Goes the Tick
  137.  
  138. In order to assess the accuracy and robustness of NTP-synchronized
  139. clocks, I use a special-purpose simulator that processes the statistical
  140. data and closely mimics the operation of the NTP daemon on a real
  141. machine. The program reads a file such as produced by the handy-dandy
  142. ASCII conversion routine described in Appendix A and produces a file
  143. showing the resulting behavior of the filtering, selecting, combining
  144. and phase-lock-loop algorithms. An example of the data input to this
  145. program is shown below (all numbers are in milliseconds, except the day:
  146. MJD and the Code: hex).
  147.  
  148.             Day      ToD  Code   Offset   Delay Disper
  149.            -------------------------------------------
  150.            48289    79369 6115          -4    39    12
  151.            48289   106624 4204          -5    63     4
  152.            48289   119282 120a      -13564   113  3623
  153.            48289   185707 120a      -13129   115  3623
  154.            48289   213995 6115          -3    38    12
  155.            48289   217210 420e          -6   116    27
  156.            48289   223560 6204          -5    62     3
  157.            48289   244875 2216          -6   192    27
  158.            48289   253305 130a        1182   153  2043
  159.            48289   276211 2208          -6   116    27
  160.            48289   277278 4209          -6   117    27
  161.            48289   279379 220b          -6   118    27
  162.            48289   315930 3113         -86   262    11
  163.            48289   318826 4206          -4    65     3
  164.            48289   319932 130a         819   154  2043
  165.            48289   348936 3115          -3    38    11
  166.            48289   358368 6204          -5    63     3
  167.            48289   386798 130a         578   153  2043
  168.            48289   454050 130a         574   154  2043
  169.            48289   484507 3115          -3    38    11
  170.            48289   493956 6204          -5    63     3
  171.            48289   589423 130a         585   153  2043
  172.            48289   619696 3115          -3    37    11
  173.            48289   629130 6204          -5    63     3
  174.            48289   763976 6204          -5    63     3
  175.            48289   774498 4207          -6    63     3
  176.            48289   787098 2216          -7   195    42
  177.            48289   858592 1113         -33   159    12
  178.            48289   859525 120a        -172   113    71
  179.            48289   889750 6115          -4    41    11
  180.            48289   898115 4204          -5    64     3
  181.            48289   911749 420f          -5   138    13
  182.            48289   912807 3114         -10    48     2
  183.            48289  1016709 120d         -28   132    16
  184.            48289  1023919 6115          -4    38    11
  185.            48289  1032084 4204          -5    62     3
  186.            48289  1056335 2216          -8   190    46
  187.  
  188. The last two hex digits of the Code indicate the peer ID, from which we
  189. learn that peer 0x0a is in awful shape, for what cause I do not know.
  190. Host 0x13 happens to be ISI fuzzball, shows a systematic offset of --33
  191. ms due to a well known asymmetric path one way via DARTnet and the other
  192. way via NSFnet. Note that these are unfiltered data; things look much
  193. better after the NTP algorithms get to work.
  194.  
  195. The following table is a snapshot taken from dcn2.udel.edu for every
  196. interface of every DARTnet router, along with fuzzball primary servers
  197. wwvb.isi.edu [128.9.2.129], dcn1.udel.edu [128.4.0.1] and dcn5.udel.edu
  198. [128.4.0.5], as well as a selected fuzzball secondary server
  199. clock.sura.net [192.80.214.42]. The data show the peer ID, IP address
  200. and (filtered) delay, offset and weight. The weight, which is used by
  201. the clock-combining algorithm (adapted from that used by NIST to produce
  202. the U.S. standard timescale from a bunch of cesium clocks), is computed
  203. as the quotient 65536 divided by the synchronization distance, itself
  204. computed as the total dispersion plus one-half the total delay to the
  205. root of the synchronization subnet, as described in the NTP Version 3
  206. specification. For this purpose higher values of weight indicate higher
  207. quality of data, as determined by the algorithm.
  208.  
  209.           ID        Address        Dly  Ofs  Wgt  
  210.           --------------------------------------
  211.           4 +  [140.173.32.2]      25   -5   1638
  212.           5 -  [140.173.16.1]      19   -6   2114
  213.           6 +  [140.173.16.2]      25   -5   1985
  214.           7 +  [140.173.64.1]      26   -4   1337
  215.           8 .  [140.173.64.2]      78   -13  840
  216.           9 .  [140.173.128.1]     78   -13  840
  217.           10 x [140.173.128.2]     80   -65  56
  218.           11 . [140.173.112.1]     78   -13  840
  219.           12 . [140.173.112.2]     88   -4   885
  220.           13 + [140.173.96.1]      87   -4   897
  221.           14 . [140.173.96.2]      77   -14  819
  222.           15 . [140.173.80.1]      90   -5   840
  223.           16 + [140.173.80.2]      87   -5   897
  224.           17 + [140.173.144.1]     88   -5   897
  225.           18 + [140.173.144.2]     90   -3   885
  226.           19 x [128.9.2.129]       169  -41  579
  227.           20 - [128.4.0.1]         39   -7   2340
  228.           21 * [128.4.0.5]         38   -3   3120
  229.           22 . [192.80.214.42]     55   -9   736
  230.  
  231. Following the peer ID is the tattletale indicator, which shows the
  232. status of that entry as determined by the NTP algorithms. The tattletale
  233. coding is as described above for the statistics file entry. The entry
  234. marked "*" is the current synchronization source as determined by NTP;
  235. however, in the case of this particular fuzzball this source is not used
  236. to synchronize the machine itself, since a primary reference source is
  237. available internally. The two entries marked "x" have been discarded by
  238. the intersection (modified Marzullo) algorithm as suspect falsetickers,
  239. while the entries marked "." have been discarded from the end of the
  240. surviving candidate list (the fuzzballs keep a maximum of ten of the
  241. lowest-distance candidates on that list). The entries marked "-" have
  242. been discarded by the clustering (Mills) algorithm, while the remaining
  243. entries marked "+" have survived all hazards and are therefore eligible
  244. for processed by the weighted-average (NIST) combining algorithm which
  245. produces the final local-clock correction (whew).
  246.  
  247. For comparison and evaluation, the following two tables show the offset,
  248. delay and dispersion for the last eight measurements made from
  249. dcn2.udel.edu to wwvb.isi.edu and dcn5.udel.edu. Note that the
  250. dispersion (due primarily to maximum assumed local-clock frequency
  251. offset) indicates the approximate age of the sample, in minutes.
  252.  
  253.      wwvb.isi.edu
  254.           Offset    -68  -41  -42  -38  -66  -69  -38  -71
  255.           Delay     227  169  174  166  224  227  164  236
  256.           Disper    1    17   34   50   67   83   117  133
  257.      
  258.      dcn5.udel.edu
  259.           Offset    -3   -4   -4   -5   -5   -6   -4   -4
  260.           Delay     38   39   38   39   39   44   38   39
  261.           Disper    1    5    9    13   17   22   26   30
  262. The systematic offset of wwvb.isi.edu is clearly apparent, as is the
  263. typical NSFnet noise on that the leg of the roundtrip path. The
  264. dcn5.udel.edu data show low dispersions typical of the quiet DARTnet.
  265. However, dcn5.udel.edu shows a systematic offset of about -5 ms due to
  266. subtle delays in the internal clock-synchronization scheme used locally.
  267. When we get all our wires uncrossed here, that curiosity should
  268. disappear. Nevertheless, what experimenters usually need is a
  269. calibration of experiment host and router offsets between each other, in
  270. order to determine a common reference timescale. In principle, it
  271. doesn't matter what the systematic offsets are, just that they can be
  272. determined retrospectively from the statistics file.
  273.  
  274. 4. How Tight the Tick
  275.  
  276. The question yet remains on the accuracy and stability of the particular
  277. local clock on each experiment host and router. It is not possible to
  278. observe this directly with the available hardware for other than
  279. fuzzballs, and, even for fuzzballs requires a precision interval counter
  280. (such as constructed by one of my students). However, the NTP local-
  281. clock model is described as an adaptive-parameter, type-II, phase-lock
  282. loop (PLL). Fortunately, one of the parameters estimated, called the
  283. compliance (Greek tau in Appendix G of the NTP Version 3 specification),
  284. is a particularly useful indicator of local-oscillator stability. Large
  285. negative values indicate an oscillator with stabilities in the range
  286. better than 1 ms per day, while more positive values indicate other
  287. oscillators varying over the range up to several seconds per day.
  288.  
  289. For fuzzballs (only) an indirect indicator of the compliance is the poll
  290. interval, which is included in all NTP packets. This value, expressed as
  291. a power of 2 in seconds, is shown in both the Unix and fuzzball
  292. monitoring displays. In version 3 of the protocol the poll interval is
  293. tied to the compliance, with effect that the interval varies from about
  294. a minute for the most unstable oscillator (the power grid) to about 17
  295. minutes for most stable oscillators (LORAN-C receivers and cesium
  296. clocks). The displayed value can vary from 6 through 10, with the
  297. highest number indicating the highest stability. Ordinarily, this tidbit
  298. would be important only for experiments designed to run for some time
  299. (hours) while NTP is turned off.
  300.  
  301. One of the goals of the NTP Version-3 models is to demonstrate
  302. accuracies to the millisecond, at least in DARTnet-like systems, as well
  303. as stabilities to the order of a millisecond per day. If these goals can
  304. be obtained in principle, a DARTnet platform could be synchronized, then
  305. NTP could be turned off during the experiment itself, if need be, with
  306. assurance the platform would not wander by more than fractional
  307. milliseconds during the experiment.
  308.  
  309. In order to test whether these goals can be met, the simulator
  310. previously described was used to process the statistics data collected
  311. from dcn2.udel.edu for the last couple of weeks. A snatch of these data
  312. covering approximately the same interval used previously to illustrate
  313. the statistics format is shown below:
  314.  
  315.        Time       Offset     Freq    Delay    Dist Clks   RefID
  316.       ---------------------------------------------------------
  317.       20001.3     -5.711   -0.0023    41.0    33.8    8      20
  318.       20003.6     -5.708   -0.0023    40.0    31.5    6      20
  319.       20005.8     -5.704   -0.0023    41.0    29.5    6      20
  320.       20008.0     -5.699   -0.0022    40.0    28.5    6      20
  321.       20010.3     -5.698   -0.0022    39.0    29.9    6      20
  322.       20012.5     -5.698   -0.0022    39.0    29.5    6      20
  323.       20014.8     -5.697   -0.0022    42.0    30.8    6      20
  324.       20017.0     -5.692   -0.0022    39.0    31.5   10      20
  325.       20019.3     -5.689   -0.0022    39.0    28.4   10      20
  326.       20021.5     -5.686   -0.0022    39.0    28.6   10      20
  327.       20023.8     -5.683   -0.0022    40.0    28.8   10      20
  328.       20026.1     -5.680   -0.0021    38.0    27.7   10      20
  329.       20028.3     -5.676   -0.0021    38.0    29.4   10      20
  330.       20030.6     -5.673   -0.0021    39.0    28.1   10      20
  331.       20032.8     -5.670   -0.0021    39.0    30.3   10      20
  332.       20035.0     -5.662   -0.0021    42.0    32.2   10      20
  333.       20037.3     -5.657   -0.0020    40.0    28.6   10      20
  334.       20039.5     -5.652   -0.0020    40.0    30.6   10      20
  335.       20041.8     -5.649   -0.0020    39.0    29.9   10      20
  336.       20044.0     -5.645   -0.0020    39.0    31.6   10      20
  337.       20046.3     -5.642   -0.0020    39.0    29.5   10      20
  338.       20048.5     -5.637   -0.0019    40.0    29.7   10      20
  339.       20050.7     -5.630   -0.0019    40.0    30.4    9      20
  340.       20053.0     -5.623   -0.0019    40.0    33.0    9      20
  341.       20055.2     -5.615   -0.0018    39.0    29.6    9      20
  342.       20057.5     -5.607   -0.0018    39.0    28.7    9      20
  343.       20059.7     -5.600   -0.0018    39.0    30.4    9      20
  344.       20062.0     -5.593   -0.0017    40.0    27.4    9      20
  345.       20064.2     -5.589   -0.0017    43.0    33.5    9      20
  346.       20066.5     -5.580   -0.0017    40.0    29.9   10      20
  347.       20068.7     -5.572   -0.0017    41.0    30.7   10      20
  348.       20070.9     -5.567   -0.0016    38.0    27.2    7      20
  349.       20075.4     -5.552   -0.0016    38.0    32.7    9      20
  350.       20084.5     -5.528   -0.0014    40.0    36.0    9      20
  351.       20102.6     -5.470   -0.0011    38.0    40.1    9      20
  352.       20120.7     -5.415   -0.0007    37.0    42.4    9      20
  353.       20138.9     -5.377   -0.0005    38.0    40.0    9      20
  354.       20159.1     -5.340   -0.0003    38.0    43.4    8      21
  355.       20161.3     -5.338   -0.0003    39.0    38.8    8      21
  356.       20163.6     -5.335   -0.0002    38.0    41.2    8      21
  357.       20165.8     -5.331   -0.0002    38.0    39.6    8      21
  358.       20168.1     -5.327   -0.0002    38.0    39.5    8      21
  359.       20170.3     -5.323   -0.0002    37.0    38.9    8      21
  360.       20174.8     -5.316   -0.0002    37.0    42.4    8      21
  361.       20177.1     -5.313   -0.0001    38.0    37.1    8      21
  362.  
  363. In this table time is shown in minutes past a certain epoch, which is
  364. not important here, while the offset, delay and distance are shown in
  365. milliseconds. The frequency column shows the computed stability in
  366. parts-per-million. The sixth column shows the number of peers remaining
  367. on the candidate list after the intersection algorithm (10 max), while
  368. the last column shows the resulting synchronization source selection.
  369.  
  370. The above data shows that, at least for periods of a couple of hours,
  371. the "real" time stays within a few hundred microseconds and the
  372. stability within a few hundred microseconds per day. However, while
  373. these data are typical of DARTnet behavior, at least for the present
  374. lightly loaded network and with the currently broken kernel clock code,
  375. performance is not always this good. Following is a summary of the
  376. network behavior over about the last couple of weeks.
  377.  
  378.       ID Samples      Mean    StdDev       Max       Min
  379.       --------------------------------------------------
  380.        0    5097    -4.122     1.719     0.983    -7.480
  381.        4     100    -4.930     7.879    73.000    -7.007
  382.        5    1608    -4.368     9.055     7.000  -103.000
  383.        6    1901    -3.967     2.810    19.000   -11.000
  384.        7    1929    -4.144     2.509    26.000   -11.000
  385.        8     558    -3.533     4.003    47.000   -26.000
  386.        9     547    -3.321     5.059    65.000   -25.000
  387.       10     968  -770.688  1025.602  1182.000-13564.000
  388.       11     549    -3.390     4.825    69.000   -25.000
  389.       12     951    -3.656     3.878    72.000   -17.000
  390.       13     445    -5.972     7.195    80.000   -30.000
  391.       14     549    -3.435     4.994    74.000   -25.000
  392.       15     945   -17.458   126.391    45.000 -1249.000
  393.       16     951    -3.561     4.118    71.000   -17.000
  394.       17     941    -3.604     3.735    72.000   -17.000
  395.       18     938    -6.643    29.347    56.000  -570.000
  396.       19    1017   -40.530    16.283    19.347  -234.000
  397.       20    2146    -5.368     4.084    27.000   -47.000
  398.       21    5419    -3.320     3.102    75.000   -24.000
  399.       22     499    -9.317    11.103    46.000   -76.000
  400.  
  401. In this table the peer ID corresponds to the fuzzball monitoring table
  402. described previously, while the statistics are computed from the offset
  403. data in the statistics file described previously. However, peer ID zero
  404. corresponds to the output of the simulation program, so represents the
  405. computed offset of the local platform (dcn2.udel.edu) relative to the
  406. mean timescale represented by all the DARTnet players. This timescale is
  407. provably the best clock in the house, with maximum error over the entire
  408. interval less than about +-3.5 ms relative to the mean. Most of the
  409. residual error is due to transients created when the machine was
  410. rebooted. Obviously, something was certainly broken at entry 10 at some
  411. time over the period of observation and probably at entry 15. The
  412. systematic offset previously mentioned with entry 19 (wwvb.isi.edu) is
  413. clearly apparent, but entry 15 doesn't look to good, either.
  414.  
  415. 5. How Timely the Tick
  416.  
  417. There is some concern about the time it takes to tune the NTP
  418. contraption exactly on frequency. The long answer is that it takes days.
  419. The short answer is the long answer stashed in a file that is read upon
  420. reboot. The xntpd program presently does this, which considerably
  421. shortens the "training" interval while the initial transient in setting
  422. the clock damps out. While the PLL is mildly underdamped and stabilizes
  423. offset within a few percent in about an hour, large offset surges can
  424. produce frequency errors that may not completely settle down until a day
  425. or two have elapsed. Ordinarily, these niceties should not be of
  426. concern, unless accuracies to the millisecond are required. Rarely do
  427. the local-clock offsets exceed ten milliseconds even right after reboot.
  428.  
  429. Since reboot frequencies may be somewhat larger than usual with DARTnet
  430. activities, some consideration of the initial synchronization procedure
  431. may help explain what may at first seem odd behavior. When an NTP client
  432. comes up it starts sending packets to each of its configured buddies.
  433. Upon receiving a validated reply (eight separate sanity checks plus
  434. crypto-authentication between fuzzballs - optional with otherballs), the
  435. offset/delay/dispersion sample is stashed in an eight-stage shift
  436. register, one for each peer.
  437.  
  438. As succeeding samples arrive, the filtering, selecting and combining
  439. engines puff away to digest the data. As the shift registers fill up the
  440. overall synchronization distance computed from the dispersion data
  441. ordinarily decreases below a preset threshold, in the case of fuzzballs
  442. 1000 ms. During all this time the algorithms are doing their best to
  443. select the best out of possibly very noisy peers and the tattletale
  444. indicators wiggle as advertised. However, the local clock is not updated
  445. until the overall distance falls below the preset threshold. At this
  446. point the PLL kicks in and the local clock begins to hum. During the
  447. initial interval after reboot one peer may tattle "#" and many (up to
  448. 10) tattle "+". Ordinarily, after about ten minutesthe bumps and grinds
  449. should eventually converge with the "#" be replaced by "*" and only a
  450. few "+".
  451.  
  452. It is important to note that, during the initial acquisition interval
  453. and others where the timing noise is exceptionally high, the number of
  454. peers apparently eligible for synchronization may seem bogusly large.
  455. The clustering algorithm used to select the best subset from a bunch of
  456. clocks compares the dispersion of each clock relative to the ensemble of
  457. all available clocks and tosses outlyers until the ensemble dispersion
  458. is less than any clock dispersion or the number of clocks hits a lower
  459. bound, currently 3. Details of this intricate procedure can be found in
  460. the NTP Version 3 specification, Section 4.2.2.
  461.  
  462. 6. How Will the Tick
  463.  
  464. There are several projects that need to be completed before the clock
  465. assemblage can be considered completely stable. First and most urgently,
  466. Van Jacobson's fixes for the SunOS timekeeping code must be installed in
  467. all routers and experiment hosts. Second, the NTP configuration at
  468. dcn1.udel.edu needs to be augmented to include the experiment hosts as
  469. the addresses are determined. Note that the watchkeeper and UDel
  470. experiment host pogo.udel.edu have routes to DARTnet and associated
  471. local nets, but not the other DCnet hosts or time servers.
  472.  
  473. A problem remains on how to synchronize the experiment hosts and
  474. possibly other hosts on the various local nets associated with DARTnet.
  475. While those hosts that know about DARTnet and its associated nets, the
  476. ISI and UDel time servers don't know about them. It is probably not a
  477. good idea to teach these servers about the DARTnet topology, since there
  478. are numerous other clients of these servers on those local nets and the
  479. preferred routes to them should not transit DARTnet. My suggested
  480. solution is to configure each of the experiment hosts to chime time
  481. exclusively from DARTnet routers. While this results the most accurate
  482. and stable time synchronization for DARTnet experiments, there may be a
  483. certain reluctance to surrender the clock to an experimental contraption
  484. in the first place. I don't have trouble doing that here and would
  485. welcome discussion on the issue.
  486.  
  487. Finally, some way must be found for the time server at ISI and the
  488. watchkeeper at UDel to schmooze with each other via DARTnet. This is
  489. critical for calibrating the ISI timescale against the UDel timescale
  490. when some DARTers chime with one and others chime with the other.
  491. Assuming a second Ethernet interface can be spliced to this beast,
  492. individual host routes can be set up so that traffic only for these
  493. hosts will ride the DARTnet rails, while traffic for other hosts on the
  494. same net will wander the NSFnet swamps. It is probably not a good idea
  495. to avoid the problem by forcing all platforms to synchronize to only one
  496. of the two servers, but not the other, since they occasionally suffer
  497. radio timewarps, are rebooted, lose power, loose connectors and whatnot.
  498. In fact, the present situation is somewhat dangerous in that only two
  499. servers with two radios are available.
  500.  
  501. The preferred way to avoid the timewarp problem is to splice a third
  502. radio to the UDel watchkeeper, so that it in effect becomes a dedicated
  503. primary server for DARTnet. I can do that if Steve Casner sends me his
  504. second WWVB clock as promised, which I can then (fix?) and install on
  505. dcn2.udel.edu. I am hoping Spectracom will give me one of their new WWVB
  506. clocks for me to evaluate, which also works. Eventually, the clock and
  507. the watchkeeper machine itself are to move to MIT; but, in that case, I
  508. assume that some way can be found to splice it to DARTnet from there.
  509.  
  510. Appendix A. Handy Programs
  511.  
  512. I do most of my data-reduction chores on a souped-up i386 with Microsoft
  513. QuickC. Among the handy-dandy programs used are two gems, one to convert
  514. the statistics file from binary format to ASCII, and the other to
  515. simulate the NTP algorithms as driven by the ASCII file. I offer them
  516. for what it's worth and not necessarily as examples of programming
  517. virtuosity.
  518.  
  519. The first program below reads an input file in the format previously
  520. described and writes its ASCII representation on the output file. It
  521. also does some error checking to cast out junk that can happen due to
  522. various kinds of server calamities, like updating software versions. The
  523. command line consists of the input file name, the output file name and a
  524. selector. If the selector matches the ID of a peer, statistics for only
  525. that peer are produced. If the selector is zero, statistics for all
  526. peers are produced.
  527.  
  528. The second program is the NTP simulator. It prompts for the name of an
  529. ASCII statistics file and then cranks the NTP algorithms to produce the
  530. data shown in the main body of this report. Output goes the control
  531. terminal or can be redirected to a file. This is the same program from
  532. which the sample C programs given in Appendix I of the NTP Version 3
  533. specification were extracted.
  534.  
  535. /* this program reads a statistics file and converts to ascii */
  536.  
  537. #include <stdio.h>
  538. #include <ctype.h>
  539. #define UTCDAY 41317            /* MJD at 1 Jan 1972 */
  540.  
  541. long int n1 = 0, n2 = 0;
  542. FILE *fopen(), *fp_in, *fp_out;
  543.  
  544. main(int argc, char *argv[]) {
  545.      int n;
  546.  
  547.      if (argc < 4) {
  548.           printf ("usage: infile outfile peer\n"); exit (1);
  549.           }
  550.      if ((fp_in = fopen (argv[1], "rb")) == NULL) {
  551.           printf ("input file does not exist\n"); exit (1);
  552.           }
  553.      fp_out = fopen (argv[2], "w");
  554.      if (sscanf(argv[3], "%i", &n) != 1) {
  555.           printf ("invalid argument\n"); exit (1);
  556.           }
  557.      convert_file (fp_in, fp_out, n);
  558.      fclose (fp_in); fclose (fp_out); return (0);
  559.      }
  560.  
  561. convert_file (FILE *input, FILE *output, int sel) {
  562.      unsigned int buf[8], i, c, date, daychk, hour, minute;
  563.      unsigned long time;
  564.      float second;
  565.      struct rec {            /* record produced by fuzzball */
  566.           unsigned int date;
  567.           unsigned long time;
  568.           unsigned int code;
  569.           long offset;
  570.           int delay;
  571.           int disp;
  572.           } *p;
  573.  
  574.      n1 = 0; n2 = 0; daychk = 0;
  575.      while (!feof(input)) {
  576.           c = fread(buf, 2, 8, input);
  577.           c = buf[1]; buf[1] = buf[2]; buf[2] = c; /* endian */
  578.           c = buf[4]; buf[4] = buf[5]; buf[5] = c;
  579.           p = buf;
  580.           date = (p->date & 0x3fff)+UTCDAY; time = p->time;
  581.           if (daychk == 0) daychk = date;
  582.           if ((date > 50000) || (date < daychk) ||
  583.                (time > 86400000) || (p->code == 0)) continue;
  584.           n1++;
  585.           if ((sel > 0) && (sel != (p->code & 0xff))) continue;
  586.           n2++;
  587.           hour = time/3600000;
  588.           minute = (time%3600000)/60000;
  589.           second = (time%60000)/1000.;
  590.      /*        fprintf(output, "%6u%3i:%02i:%06.3f", date,
  591.                hour, minute, second);                  */
  592.           fprintf(output, "%6u%9lu", date, time);
  593.           fprintf(output, " %04x%12li%6i%6i\n", p->code, p->offset,
  594.                p->delay, p-> disp);
  595.           }
  596.      printf("input %li  output %li\n", n1, n2);
  597.      }
  598.  
  599. /*
  600.    NTP simulator
  601.  
  602.    This program reads sample data from a history file and simulates the
  603.    behavior of the NTP filter, selection and local-clock algorithms.
  604.  */
  605. #include <stdio.h>
  606. #include <ctype.h>
  607. #include <math.h>
  608. /*
  609.    Parameters
  610.  */
  611. #define C_SYNC 0x8000         /* synch okay */
  612. #define NMAX 40               /* max clocks */
  613. #define FMAX 8                /* max filter size */
  614. #define DAY 86400.            /* one day of seconds */
  615. #define MAXREACH 8192.        /* reachability timeout */
  616. #define HZ 1000.              /* clock rate */
  617. #define FILTER .5             /* filter weight */
  618. #define SELECT .75            /* select weight */
  619. #define MAXSTRAT 15.          /* max stratum */
  620. #define MAXSKEW 1.            /* max skew error per MAXAGE */
  621. #define MAXAGE 86400.         /* max clock age */
  622. #define MAXDISP 16.           /* max dispersion */
  623. #define MAXCLOCK 10           /* max candidate clocks */
  624. #define MINCLOCK 3            /* min survivor clocks */
  625. /*
  626.    Function declarations
  627.  */
  628. void filter(int, double, float, float); /* filter algorithm */
  629. int dts(void);                /* intersection algorithm */
  630. int select(void);             /* selection algorithm */
  631. float combine(void);          /* combining algorithm */
  632. float distance(int);          /* compute synch distance */
  633. void stat(int, float);        /* update statistics */
  634. void display(void);           /* display statistics */
  635. void reset(void);             /* reset statistics */
  636. void exit(int);               /* off the bus */
  637. /*
  638.    Peer state variables
  639.  */
  640. float filtp[NMAX][FMAX];      /* filter offset samples */
  641. float fildp[NMAX][FMAX];      /* filter delay samples */
  642. float filep[NMAX][FMAX];      /* filter dispersion samples */
  643. float tp[NMAX];               /* offset */
  644. float dp[NMAX];               /* delay */
  645. float ep[NMAX];               /* dispersion */
  646. float rp[NMAX];               /* last offset */
  647. double utc[NMAX];             /* last timestamp */
  648. int cd[NMAX];                 /* code */
  649. int st[NMAX];                 /* stratum */
  650. /*
  651.    System state constants/variables
  652.  */
  653. float rho = 1./HZ;            /* max reading error */
  654. float phi = MAXSKEW/MAXAGE;   /* max skew rate  */
  655. float maxdist = MAXDISP/16.;  /* max distance */
  656. float bot, top;               /* confidence interval limits */
  657. float theta;                  /* clock offset */
  658. float delta;                  /* roundtrip delay */
  659. float epsil;                  /* dispersion */
  660. int nclock;                   /* survivor clocks */
  661. int source;                   /* clock source */
  662. int n1, n2;                   /* min/max clock ids */
  663. /*
  664.    Local clock constants
  665.  */
  666. double Kf = 4194394.;         /* frequency weight */
  667. double Kg = 256.;             /* phase weight */
  668. double Kh = 8192.;            /* compliance weight */
  669. double Ks = 16.;              /* compliance maximum */
  670. double Kt = 8192.;            /* compliance multiplier */
  671. double sigma = 4.;            /* adjustment interval */
  672. /*
  673.    Local clock variables
  674.  */
  675. double Vs;                    /* clock filter output */
  676. double tau;                   /* PLL time constant */
  677. double mu;                    /* update interval */
  678. double freq;                  /* frequency error */
  679. double phase;                 /* phase error */
  680. double comp;                  /* error estimate */
  681. double tstamp;                /* at the tone */
  682. /*
  683.    Statistics vector
  684.  */
  685. double d0[NMAX];              /* sample count */
  686. double d1[NMAX];              /* sum of samples */
  687. double d2[NMAX];              /* sum of squares of samples */
  688. double d3[NMAX];              /* sample maximum */
  689. double d4[NMAX];              /* sample minimum */
  690. /*
  691.    Input file format
  692.  */
  693. unsigned int date;            /* day number (MJD) */
  694. unsigned long timex;          /* time of day (ms past midnight) */
  695. unsigned int code;            /* entry code */
  696. float offset;                 /* clock offset (ms) */
  697. float delay;                  /* total delay (ms) */
  698. float disp;                   /* total dispersion (ms) */
  699. /*
  700.    Temporary lists
  701.  */
  702. float list[3*NMAX];           /* temporary list */
  703. int index[3*NMAX];            /* index list */
  704. /*
  705.    Stuff
  706.  */
  707. unsigned int datey;           /* starting date */
  708. FILE *fopen(), *fp_in;        /* input file handle */
  709. char fn_in[50];               /* input file name */
  710. float t, u, v, w, x;          /* float temps */
  711. int i, j, k, m, n;            /* int temps */
  712.  
  713. main() {
  714.      printf ("input file: "); scanf ("%s", fn_in);
  715.      fp_in = fopen (fn_in, "r");
  716.      if (fp_in == NULL) {
  717.           printf ("input file %s does not exist\n", fn_in);
  718.           exit (1);
  719.           }
  720.      n1 = 1; n2 = 25;
  721.      reset();
  722.      printf("    Time     Offset      Freq   Delay    Disp    Dist  Clk
  723. Src\n");
  724. /*
  725.    Read next sample
  726.  */
  727.      while (fscanf(fp_in, "%i%lu%x%f%f%f",
  728.           &date, &timex, &code, &offset, &delay, &disp) != EOF) {
  729.           timex = timex/1e3; offset = offset/1e3;
  730.           delay = delay/1e3; disp = disp/1e3;
  731.           if (datey <= 0) datey = date;
  732.           tstamp = (date-datey)*DAY+timex;
  733.           if (utc[0] <= 0) utc[0] = tstamp;
  734.           for (n = n1; n <= n2; n++) {
  735.                if (utc[n] <= 0) utc[n] = tstamp;
  736.                if ((cd[n] != 0) && ((tstamp-utc[n]) > MAXREACH)) {
  737.                     for (i = FMAX-1; i >= 0; i--) {
  738.                          filtp[n][i] = 0.; fildp[n][i] = 0.;
  739.                          filep[n][i] = MAXDISP;
  740.                          }
  741.                     tp[n] = 0.; rp[n] = 0.;
  742.                     ep[n] = MAXDISP; utc[n] = 0.;
  743.                     cd[n] = 0;
  744.                     }
  745.                }
  746.           n = code & 0xff;
  747.           if ((n < n1) || (n > n2)) continue;
  748.           cd[n] = 0;
  749.           if (code&C_SYNC == 0) continue;
  750.           cd[n] = code; st[n] = (code&0x0f00)/256.;
  751. /*
  752.    Update clock filter and statistics
  753.  
  754.    n = peer id, offset = clock offset, delay = roundtrip delay to root,
  755.    disp = dispersion to root
  756.  */
  757.           t = tstamp-utc[n]; epsil = rho+phi*t+disp;
  758.           filter(n, offset, delay, epsil);
  759.           if (t > 0.) rp[n] = (tp[n]-u)/t;
  760.           stat(n, tp[n]*1e3);
  761.  /*
  762.    Determine survivors and combine
  763.  
  764.    nclock = number of survivor clocks
  765.  */
  766.           i = dts(); if (i < nclock/2) continue;
  767.           i = select(); if (i < 1) continue;
  768.           nclock = i; theta = combine();
  769.           if ((source != n) || (epsil+delta/2. >= maxdist))
  770.                continue;
  771.  /*
  772.    Update local clock and statistics
  773.  */
  774.           dp[0] = delta; ep[0] = epsil; Vs = theta;
  775.           mu = tstamp-utc[0]; utc[0] = tstamp; u = tp[0];
  776.           x = mu; if (x > 1024.) x = 1024.;
  777.           freq = freq+Vs*x/(tau*tau); phase = Vs/tau;
  778.           tau = Ks-fabs(comp); if (tau < 1.) tau = 1.;
  779.           comp = comp+(Kt*Vs-comp)*tau/Kh;
  780.           i = mu/sigma; x = 1.-1./Kg; x = 1.-pow(x, i);
  781.           tp[0] = tp[0]+x*phase+i/Kf*freq;
  782.  
  783.           printf("%8.1f%11.3f%10.4f%8.1f%8.1f%8.1f%5i%5i\n",
  784.                tstamp/60., tp[0]*1e3, 1./(sigma*Kf)*freq*1e6,
  785.                dp[0]*1e3, ep[0]*1e3, distance(0)*1e3, nclock,
  786.                cd[source]&0xff);
  787.           if (mu > 0.) rp[0] = (tp[0]-u)/mu;
  788.           stat(0, tp[0]*1e3);
  789.           }
  790.      fclose(fp_in); printf("Summary\n"); display(); exit(0);
  791.      }
  792. /*
  793.    Subroutines
  794.  
  795.    Clock filter algorithm
  796.  
  797.    n = peer id, offset = sample offset, delay = sample delay,
  798.    disp = sample dispersion; computes tp[n] = peer offset,
  799.    dp[n] = peer delay, ep[n] = peer dispersion
  800.  */
  801. void filter(int n, double offset, float delay, float disp) {
  802.      int i, j, k, m;          /* int temps */
  803.      float x;                 /* float temps */
  804.  
  805.      for (i = FMAX-1; i > 0; i--) {     /* update/shift filter */
  806.           filtp[n][i] = filtp[n][i-1];
  807.           fildp[n][i] = fildp[n][i-1];
  808.           filep[n][i] = filep[n][i-1]+phi*(tstamp-utc[n]);
  809.           }
  810.      utc[n] = tstamp; filtp[n][0] = offset-tp[0];
  811.      fildp[n][0] = delay; filep[n][0] = disp;
  812.      m = 0;                        /* construct/sort temp list */
  813.      for (i = 0; i < FMAX; i++) {
  814.           if (filep[n][i] >= MAXDISP) continue;
  815.           list[m] = filep[n][i]+fildp[n][i]/2.; index[m] = i;
  816.           for (j = 0; j < m; j++) {
  817.                if (list[j] > list[m]) {
  818.                     x = list[j]; k = index[j];
  819.                     list[j] = list[m]; index[j] = index[m];
  820.                     list[m] = x; index[m] = k;
  821.                     }
  822.                }
  823.           m = m+1;
  824.           }
  825.      if (m <= 0)              /* compute filter dispersion */
  826.           ep[n] = MAXDISP;
  827.      else {
  828.           ep[n] = 0;
  829.           for (i = FMAX-1; i >= 0; i--) {
  830.                if (i < m) {
  831.                     x = fabs(filtp[n][index[0]]-
  832.                          filtp[n][index[i]]);
  833.                     }
  834.                else x = MAXDISP;
  835.                ep[n] = FILTER*(ep[n]+x);
  836.                }
  837.           i = index[0]; ep[n] = ep[n]+filep[n][i];
  838.           tp[n] = filtp[n][i]; dp[n] = fildp[n][i];
  839.           }
  840.      return;
  841.      }
  842. /*
  843.    Intersection algorithm
  844.  
  845.    computes nclock = candidate clocks, bot = lowpoint, top = highpoint;
  846.    returns number of survivor clocks
  847. */
  848. int dts() {
  849.      int f;                   /* intersection ceiling */
  850.      int end;                 /* endpoint counter */
  851.      int clk;                 /* falseticker counter */
  852.      int i, j, k, n;          /* int temps */
  853.      float x, y;              /* float temps */
  854.  
  855.      nclock = 0; i = 0;
  856.      for (n = n1; n <= n2; n++) {  /* construct endpoint list */
  857.           if (ep[n] >= MAXDISP) continue;
  858.           nclock++;
  859.           list[i] = tp[n]-distance(n); index[i] = -1;  /* lowpoint */
  860.           for (j = 0; j < i; j++) {
  861.                if ((list[j] > list[i]) ||
  862.                     ((list[j] == list[i]) &&
  863.                     (index[j] > index[i]))) {
  864.                     x = list[j]; k = index[j];
  865.                     list[j] = list[i]; index[j] = index[i];
  866.                     list[i] = x; index[i] = k;
  867.                     }
  868.                }
  869.           i = i+1;
  870.           list[i] = tp[n]; index[i] = 0;               /* midpoint */
  871.           for (j = 0; j < i; j++) {
  872.                if ((list[j] > list[i]) ||
  873.                     ((list[j] == list[i]) &&
  874.                     (index[j] > index[i]))) {
  875.                     x = list[j]; k = index[j];
  876.                     list[j] = list[i]; index[j] = index[i];
  877.                     list[i] = x; index[i] = k;
  878.                     }
  879.                }
  880.           i = i+1;
  881.           list[i] = tp[n]+distance(n); index[i] = 1;   /* highpoint */
  882.           for (j = 0; j < i; j++) {
  883.                if ((list[j] > list[i]) ||
  884.                     ((list[j] == list[i]) &&
  885.                     (index[j] > index[i]))) {
  886.                     x = list[j]; k = index[j];
  887.                     list[j] = list[i]; index[j] = index[i];
  888.                     list[i] = x; index[i] = k;
  889.                     }
  890.                }
  891.           i = i+1;
  892.           }
  893.      if (nclock <= 0) return 0;
  894.      for (f = 0; ; f++) {     /* find low endpoint */
  895.           clk = 0; end = 0;
  896.           for (j = 0; j < i; j++) {
  897.                end = end-index[j]; bot = list[j];
  898.                if (end >= (nclock-f)) break;
  899.                if (index[j] == 0) clk = clk+1;
  900.                }
  901.           end = 0;
  902.           for (j = i-1; j >= 0; j--) {
  903.                end = end+index[j]; top = list[j];
  904.                if (end >= (nclock-f)) break;
  905.                if (index[j] == 0) clk = clk+1;
  906.                }
  907.           if (clk <= f) break;
  908.           }
  909.      return nclock-clk;
  910.      }
  911. /*
  912.    Selection algorithm
  913.  
  914.    bot = lowpoint, top = highpoint; computes nclock = number of
  915.    candidate clocks, index = candidate index list, source = clock
  916.    source, theta = clock offset, delta = roundtrip delay,
  917.    epsil = dispersion; returns number of survivor clocks
  918. */
  919. int select() {
  920.      float xi;                /* max select dispersion */
  921.      float eps;               /* min peer dispersion */
  922.      int i, j, k, m, n;       /* int temps */
  923.      float x, y, z;           /* float temps */
  924.  
  925.      m = 0;
  926.      for (n = n1; n <= n2; n++) {  /* construct/sort candidate list */
  927.           if ((st[n] > 0) && (st[n] < MAXSTRAT) &&
  928.                (tp[n] >= bot) && (tp[n] <= top)) {
  929.                list[m] = MAXDISP*st[n]+distance(n); index[m] = n;
  930.                for (j = 0; j < m; j++) {
  931.                     if (list[j] > list[m]) {
  932.                          x = list[j]; k = index[j];
  933.                          list[j] = list[m];
  934.                          index[j] = index[m];
  935.                          list[m] = x; index[m] = k;
  936.                          }
  937.                     }
  938.                m = m+1;
  939.                }
  940.           }
  941.      nclock = m;
  942.      if (m == 0) {
  943.           source = 0; return m;
  944.           }
  945.      if (m > MAXCLOCK) m = MAXCLOCK;
  946.      while (1) {                   /* cast out falsetickers */
  947.           xi = 0.; eps = MAXDISP;
  948.           for (j = 0; j < m; j++) {
  949.                x = 0.;
  950.                for (k = m-1; k >= 0; k--)
  951.                     x = SELECT*(x+fabs(tp[index[j]]-
  952.                          tp[index[k]]));
  953.                if (x > xi) {       /* max(xi) */
  954.                     xi = x; i = j;
  955.                     }
  956.                x = ep[index[j]]+phi*(tstamp-utc[index[j]]);
  957.                if (x < eps) eps = x;   /* min(eps) */
  958.                }
  959.           if ((xi <= eps) || (m <= MINCLOCK)) break;
  960.           if (index[i] == source) source = 0;
  961.           for (j = i; j < m-1; j++) index[j] = index[j+1];
  962.           m = m-1;
  963.           }
  964.      i = index[0];                 /* declare winner */
  965.      if (source != i)
  966.           if (source == 0) source = i;
  967.           else if (st[i] < st[source]) source = i;
  968.      theta = tp[i]; delta = dp[i];
  969.      epsil = ep[i]+phi*(tstamp-utc[i])+xi+fabs(tp[i]);
  970.      return m;
  971.      }
  972. /*
  973.    Combining algorithm
  974.  
  975.    index = candidate index list, m = number of candidates;
  976.    returns combined clock offset
  977.  */
  978. float combine() {
  979.      int i;                   /* int temps */
  980.      float x, y, z;           /* float temps */
  981.  
  982.      z = 0. ; y = 0.;
  983.      for (i = 0; i < nclock; i++) {     /* compute weighted offset */
  984.           j = index[i]; x = distance(j);
  985.           z = z+tp[j]/x; y = y+1./x;
  986.           }
  987.      return z/y;              /* normalize */
  988.      }
  989. /*
  990.    Compute synch distance
  991.  
  992.    n = peer id; returns synchronization distance
  993.  */
  994. float distance(int n) {
  995.  
  996.      return ep[n]+phi*(tstamp-utc[n])+dp[n]/2.;
  997.      }
  998. /*
  999.    Update statistics
  1000.  */
  1001. void stat(int n, float u) {
  1002.  
  1003.      d0[n] = d0[n]+1; d1[n] = d1[n]+u; d2[n] = d2[n]+u*u;
  1004.      if (u > d3[n]) d3[n] = u; if (u < d4[n]) d4[n] = u;
  1005.      return;
  1006.      }
  1007. /*
  1008.    Display statistics
  1009.  */
  1010. void display() {
  1011.      int i;                   /* int temps */
  1012.  
  1013.      printf(" ID    Samp      Mean    StdDev       Max       Min\n");
  1014.      for (i = 0; i < NMAX; i++) {
  1015.           if (d0[i] <= 1.) continue;
  1016.           v = sqrt(d2[i]/d0[i]-d1[i]/d0[i]*d1[i]/d0[i]);
  1017.           /*   v = sqrt(d2[i]/(2*(d0[i]-1.)));   */
  1018.           printf("%3i%8.0f%10.3f%10.3f%10.3f%10.3f\n",
  1019.                i, d0[i], d1[i]/d0[i], v, d3[i], d4[i]);
  1020.           /*   fprintf(fp_out,
  1021.                     "%3i%8.0f%10.3f%10.3f%12.3f%12.3f%10.3f\n",
  1022.                     i, d0[i], d1[i]/d0[i], v, d3[i], d4[i]);      */
  1023.           }
  1024.      }
  1025. /*
  1026.    Reset statistics
  1027.  */
  1028. void reset() {
  1029.      int i, j;                /* int temps */
  1030.  
  1031.      freq = 0.; phase = 0.; comp = 0.; tau = 1.; source = 0;
  1032.      for (i = 0; i < NMAX; i++) {
  1033.           for (j = FMAX-1; j >= 0; j--) {
  1034.                filtp[i][j] = 0.; fildp[i][j] = 0.;
  1035.                filep[i][j] = MAXDISP;
  1036.                }
  1037.           tp[i] = 0.; rp[i] = 0.; ep[i] = MAXDISP;
  1038.           utc[i] = 0.; cd[i] = 0; st[i] = 0;
  1039.           d0[i] = 0.; d1[i] = 0.; d2[i] = 0.;
  1040.           d3[i] = -1e10; d4[i] = 1e10;
  1041.           }
  1042.      }
  1043.